C
C     ..................................................................
C
C     SUBROUTINE DDBAR
C
C     PURPOSE
C	 TO COMPUTE, AT A GIVEN POINT X, AN APPROXIMATION Z TO THE
C	 DERIVATIVE OF AN ANALYTICALLY GIVEN FUNCTION FCT THAT IS 11-
C	 TIMES DIFFERENTIABLE IN A DOMAIN CONTAINING A CLOSED INTERVAL -
C	 THE SET OF T BETWEEN X AND X+H (H POSITIVE OR NEGATIVE) - USING
C	 FUNCTION VALUES ONLY ON THAT INTERVAL.
C
C      USAGE
C	 CALL DDBAR(X,H,IH,FCT,Z,)
C	 PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT
C
C     DESCRIPTION OF PARAMETERS
C	 X   - THE POINT AT WHICH THE DERIVATIVE IS TO BE COMPUTED
C	       X IS IN DOUBLE PRECISION
C	 H   - THE NUMBER THAT DEFINES THE CLOSED INTERVAL WHOSE END-
C	       POINTS ARE X AND X+H (SEE PURPOSE)
C	       H IS IN SINGLE PRECISION
C	 IH  - INPUT PARAMETER (SEE REMARKS AND METHOD)
C	       IH NON-ZERO - THE SUBROUTINE GENERATES THE INTERNAL
C			     VALUE HH
C	       IH    =	 0 - THE INTERNAL VALUE HH IS SET TO H
C	 FCT - THE NAME OF THE EXTERNAL DOUBLE PRECISION FUNCTION
C	       SUBPROGRAM THAT WILL GENERATE THE NECESSARY FUNCTION
C	       VALUES.
C	 Z   - RESULTING DERIVATIVE VALUE - DOUBLE PRECISION
C
C     REMARKS
C	 (1)  IF H = 0, THEN THERE IS NO COMPUTATION.
C	 (2)  THE (MAGNITUDE OF THE) INTERNAL VALUE HH, WHICH IS DETER-
C	      MINED ACCORDING TO IH, IS THE MAXIMUM STEP-SIZE USED IN
C	      THE COMPUTATION OF THE ONE-SIDED DIVIDED DIFFERENCES (SEE
C	      METHOD.)	IF IH IS NON-ZERO, THEN THE SUBROUTINE GENERATES
C	      HH ACCORDING TO CRITERIA THAT BALANCE ROUND-OFF AND TRUN-
C	      CATION ERROR.  HH ALWAYS HAS THE SAME SIGN AS H AND IT IS
C	      ALWAYS LESS THAN OR EQUAL TO THE MAGNITUDE OF H IN AB-
C	      SOLUTE VALUE, SO THAT ALL COMPUTATION OCCURS IN THE CLOSED
C	      INTERVAL DETERMINED BY H.
C
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C	 THE EXTERNAL FUNCTION SUBPROGRAM FCT(T) MUST BE FURNISHED BY
C	 THE USER. FCT(T) IS IN DOUBLE PRECISION
C
C     METHOD
C	 THE COMPUTATION OF Z IS BASED ON RICHARDSON'S AND ROMBERG'S
C	 EXTRAPOLATION METHOD AS APPLIED TO THE SEQUENCE OF ONE-SIDED
C	 DIVIDED DIFFERENCES ASSOCIATED WITH THE POINT PAIRS
C	 (X,X+(K*HH)/10)K=1,...,10.  (SEE FILLIPI, S. AND ENGELS, H.,
C	 ALTES UND NEUES ZUR NUMERISCHEN DIFFERENTIATION, ELECTRONISCHE
C	 DATENVERARBEITUNG, ISS. 2 (1966), PP. 57-65.)
C
C     ..................................................................
C
      SUBROUTINE DDBAR(X,H,IH,FCT,Z)
C
C
      DIMENSION AUX(10)
      DOUBLE PRECISION X,FCT,Z,AUX,A,B,C,D,DH,HH
C
C	 NO ACTION IN CASE OF ZERO INTERVAL LENGTH
      IF(H)1,17,1
C
C	 GENERATE STEPSIZE HH FOR DIVIDED DIFFERENCES
    1 C=ABS(H)
      B=H
      D=X
      D=FCT(D)
      IF(IH)2,9,2
    2 HH=.5D-2
      IF(C-HH)3,4,4
    3 HH=B
    4 HH=DSIGN(HH,B)
      Z=DABS((FCT(X+HH)-D)/HH)
      A=DABS(D)
      HH=1.D-2
      IF(A-1.D0)6,6,5
    5 HH=HH*A
    6 IF(Z-1.D0)8,8,7
    7 HH=HH/Z
    8 IF(HH-C)10,10,9
    9 HH=B
   10 HH=DSIGN(HH,B)
C
C	 INITIALIZE DIFFERENTIATION LOOP
      Z=(FCT(X+HH)-D)/HH
      J=10
      JJ=J-1
      AUX(J)=Z
      DH=HH/DFLOAT(J)
      DZ=1.E38
C
C	 START DIFFERENTIATION LOOP
   11 J=J-1
      C=J
      HH=C*DH
      AUX(J)=(FCT(X+HH)-D)/HH
C
C	 INITIALIZE EXTRAPOLATION LOOP
      D2=1.E38
      B=0.D0
      A=1.D0/C
C
C        START EXTRAPOLATION LOOP
      DO 12 I=J,JJ
      D1=D2
      B=B+A
      HH=(AUX(I)-AUX(I+1))/B
      AUX(I+1)=AUX(I)+HH
C
C        TEST ON OSCILLATING INCREMENTS
      D2=DABS(HH)
      IF(D2-D1)12,13,13
   12 CONTINUE
C        END OF EXTRAPOLATION LOOP
C
C        UPDATE RESULT VALUE Z
      I=JJ+1
      GO TO 14
   13 D2=D1
      JJ=I
   14 IF(D2-DZ)15,16,16
   15 DZ=D2
      Z=AUX(I)
   16 IF(J-1)17,17,11
C        END OF DIFFERENTIATION LOOP
C
   17 RETURN
      END
UE
C        END OF EXTRAPOLATION LOOP
C
C        UPDATE RESULT VALUE Z
      I=JJ+1
